
% path
addpath('sym') ;
addpath('ss') ;
addpath('../code') ;


% steady state
% ssfold = 'ss';
% disp('Steady-state computation...')
% tic;
% cd ss
getssnew
transfer
% cd ..
% disp(toc);

% linearized model
[F,G,Fz,Fzp,Fnp,Gz,Gzp,Gnp, ...
    GFzz,GFzpz,GFnpnp,GFznp,GFzpnp,GFzpzp, ...
    GFzpzpzp,GFzpzpnp,GFzpzpz,GFzpzz,GFzpznp,GFzpnpnp,GFzzz,GFzznp,GFznpnp,GFnpnpnp, ...
    ns,z_vec,zp_vec,np_vec] = model('order',3) ;

nv

Onn = 0.0072^2 ;

ss = 3 ; % shock size in standard deviations
time = 2:26 ;

% first order solution:
fprintf(1,'Testing first order solution...\n') ;
[PI,PHIs,PHIn] = jkl(nGzp,nFzp,nGz,nFz,nGnp,nFnp,Onn) ;

ir1 = irf(PI,PHIs,PHIn,'shock=',ss,'T=',41) ;
fprintf(1,'\n') ;

% second order solution:
fprintf(1,'Testing second order solution...\n') ;
[PI,PHIs,PHIn, ...
    Ed2z,T2zq,T2z,T2qq,T2qv,T2q] ...
    = jkl(nGzp,nFzp,nGz,nFz,nGnp,nFnp, ...
    nGFzz,nGFzpz,nGFnpnp,nGFznp,nGFzpnp,nGFzpzp, ...
    Onn) ;

ir2 = irf(PI,PHIs,PHIn,...
    T2zq,T2z,T2qq,T2qv,T2q,'shock=',ss,'T=',41) ;

fprintf(1,'\n') ;

% third order solution:
fprintf(1,'Testing third order solution...\n') ;
[PI,PHIs,PHIn, ...
    Ed2z,T2zq,T2z,T2qq,T2qv,T2q, ...
    Ed3z,T3zq,T3z,T3qq,T3qv,T3q] ...
    = jkl(nGzp,nFzp,nGz,nFz,nGnp,nFnp, ...
    nGFzz,nGFzpz,nGFnpnp,nGFznp,nGFzpnp,nGFzpzp, ...
    nGFzpzpzp,nGFzpzpz,nGFzpzpnp,nGFzpzz,nGFzpznp,nGFzpnpnp,nGFzzz,nGFzznp,nGFznpnp,nGFnpnpnp,...
    Onn) ;

ir3 = irf(PI,PHIs,PHIn,...
    T2zq,T2z,T2qq,T2qv,T2q,...
    T3zq,T3z,T3qq,T3qv,T3q,'shock=',ss,'T=',41) ;
fprintf(1,'\n') ;

varsToPlot = [1 3 14 15] ;
figure(1) ;
clf ;
for ii = 1:4
    subplot(2,2,ii) ;
    plot(time-2, 100*[ir1(varsToPlot(ii),time)/eval(z_vec(varsToPlot(ii)));...
        ir2(varsToPlot(ii),time)/eval(z_vec(varsToPlot(ii)));...
        ir3(varsToPlot(ii),time)/eval(z_vec(varsToPlot(ii)))]')
    xlabel('t') ;
    ylabel(char(z_vec(varsToPlot(ii)))) ;
    legend('o1','o2','o3') ;
end

% multi period irfs
fprintf(1,'Testing multi period shock functionality...\n') ;

sms = [ss ss] ;

ir1m = irf(PI,PHIs,PHIn,'shock=',sms,'T=',41) ;

ir2m = irf(PI,PHIs,PHIn,...
    T2zq,T2z,T2qq,T2qv,T2q,'shock=',sms,'T=',41) ;

ir3m = irf(PI,PHIs,PHIn,...
    T2zq,T2z,T2qq,T2qv,T2q,...
    T3zq,T3z,T3qq,T3qv,T3q,'shock=',sms,'T=',41) ;
fprintf(1,'\n') ;

figure(2) ;
clf ;
for ii = 1:4
    subplot(2,2,ii) ;
    plot(time-2, 100*[ir1(varsToPlot(ii),time)/eval(z_vec(varsToPlot(ii)));...
        ir2(varsToPlot(ii),time)/eval(z_vec(varsToPlot(ii)));...
        ir3(varsToPlot(ii),time)/eval(z_vec(varsToPlot(ii)))]')
    xlabel('t') ;
    ylabel(char(z_vec(varsToPlot(ii)))) ;
    legend('o1','o2','o3') ;
end

